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ABSTRACT 


This  study  is  aimed  at  development  and  application  of  new  wave  propagation  and  modeling  methods  for 
regional  waves  in  heterogeneous  crustal  waveguides  using  one-way  wave  approximation.  The  half-space 
generalized  screen  propagators  (GSP)  for  both  SH  and  P-SV  waves  have  taken  the  free  surface  into  the 
formulation  and  adopt  a  fast  dual  domain  implementation.  The  method  is  several  orders  of  magnitude  faster 
than  finite -difference  method  with  a  similar  accuracy  for  certain  problems.  It  has  been  used  for  the 
simulation  of  wave  propagation  for  high-frequency  waves  (l-25Hz)  to  a  regional  distance  (greater  than 
1000km). 

In  this  year,  we  further  develop  the  method  in  two  fronts.  First,  we  extend  the  SH  GSP  method  to  treat 
irregular  surface  topography  by  incorporating  a  coordinate  transform  into  the  method.  It  is  demonstrated 
that  our  new  approach  has  superior  wide-angle  response  to  surface  topography  over  the  PE  (parabolic 
equation)  method.  The  efficiency  of  the  screen  propagator  approach  makes  it  very  promising  for  long 
distance  Lg  simulation.  For  a  test  model  with  a  propagation  distance  of  250km  and  a  dominant  frequency 
of  1Hz,  the  screen  method  took  about  35  minutes,  while  the  boundary  element  method  took  about  72  hours. 
For  longer  propagation  distances  and  higher  frequencies,  the  factor  of  saving  could  be  huge. 

Second,  we  developed  a  P-SV  screen  propagator  for  crustal  waveguides  with  flat  surfaces.  This  is  the  major 
task  of  this  year.  Free  surface  reflection  and  conversion  have  been  incorporated  into  the  screen  propagator 
theory.  Numerical  tests  have  been  done  against  the  wavenumber  integration  and  finite -difference 
calculations.  The  results  demonstrated  the  feasibility  of  the  approach. 

Finally,  numerical  simulations  have  been  conducted  for  various  crustal  waveguide  structures,  including 
deterministic  structures,  small-scale  random  heterogeneities  and  random  rough  surfaces.  Influences  of 
random  heterogeneities  and  rough  surfaces  on  Lg  amplitude  attenuation  and  Lg  coda  formation  have  been 
studied. 


INSTRUCTION 


Short-period  regional  phases  ( Pn ,  Pg  ,  Sn ,  Lg ,  etc.)  play  an  important  role  for  both  discrimination  and 

yield  estimation  procedures  for  monitoring  the  Comprehensive  Nuclear-Test-Ban  Treaty  (CTBT). 
Amplitude  ratios  of  short-period  P  wave  energy  ( Pn  ,  P  )  to  S  wave  energy  (  Sn ,  Lg )  have  emerged  as 

particularly  promising  discriminants,  especially  at  frequencies  above  3Hz  (Walter  et  al.,  1995;  Taylor, 
1996;  Hartse  et  al.,  1997;  Fan  and  Lay,  1998).  Regional  phases  mainly  propagate  in  the  crustal  waveguide 
and  the  uppermost  mantle.  A  wide  variety  of  observational  and  theoretical  studies  show  that  the  structure 
and  multi-scale  heterogeneities  of  the  crustal  wave  guide  and  uppermost  mantle  strongly  affect  the  regional 
phases  in  waveform  properties  and  amplitudes  (Kennett,  1986,  1989;  Zhang  et  al.,  1994;  Fan  and  Lay, 
1998).  However,  due  to  the  complexity  involved  in  regional  phase  propagation  and  limit  of  capability  of 
the  existing  analysis  and  simulation  approaches,  the  mechanism  of  excitation  and  propagation  of  regional 
phases  is  still  not  appreciated.  In  these  investigations,  numerical  approaches  for  simulating  regional  phases 
are  very  useful. 
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In  the  crustal  waveguide  environment,  major  wave  energy  is  carried  by  forward  propagating  waves, 
including  forward  scattered  waves,  therefore  neglecting  backscattered  waves  in  the  propagation  will  not 
change  the  main  features  of  regional  waves  in  most  cases.  Based  on  the  concept,  a  generalized  screen 
propagator  method  (GSP)  based  on  one-way  wave  equation  has  been  developed  by  Wu  (1994),  Wu,  Jin  and 
Xie  (1999a)  and  has  been  successfully  used  to  simulate  SH  Lg  waves  in  the  complex  crustal  waveguide  and 
investigate  the  energy  partitioning  of  Lg  waves  (Wu,  Jin  and  Xie,  1999b).  Wu,  Xie  and  Wu  (1999)  further 
extended  the  GSP  method  to  treat  complex  crustal  waveguides  with  irregular  or  rough  topography  by 
incorporating  a  surface  topographic  transformation  into  the  method.  For  a  test  model  with  propagation 
distance  of  250km,  dominant  operating  frequency  of  1Hz,  the  BE  (boundary  element)  method  took  about 
72  hours,  while  the  screen  method  took  only  35  minutes  with  the  same  accuracy.  For  longer  propagation 
distances  and  higher  frequencies,  the  factor  of  saving  could  be  huge. 

Based  on  the  success  of  Lg  propagation  using  SH  screen  propagators,  in  this  paper  we  extend  the  screen 
propagator  method  to  P-SV  Lg  wave  propagation.  As  the  first  step  we  apply  the  complex  screen  method  to 
a  flat,  heterogeneous  half  space.  Unlike  SH  wave,  the  image  method  can  not  be  used  to  account  for  the 
effect  of  free  surface.  We  incorporate  the  reflection  coefficient  calculations  of  elastic  waves  on  free  surface 
into  the  complex  screen  method  for  the  effect  of  free  surface  on  P-SV  waves.  Comparison  of  synthetic 
seismograms  calculated  by  the  complex  screen  method  and  wavenumber  integration  method  for  an  elastic 
half-space  shows  feasibility  and  validity  of  the  approach.  Synthetic  seismograms  for  elastic  layered  model 
and  a  laterally  varying  crustal  waveguide  are  conducted  and  compared  with  wavenumber  integration 
method  and  finite  difference  method,  respectively.  The  results  demonstrate  that  the  complex  screen  method 
can  be  used  to  model  regional  phases  containing  Pn ,  Pg  ,  Sn  and  Lg  . 

METHODOLOGY 


Figure  1  illustrates  the  principle  of  the  method.  The  screen  method  is  a  forward  matching  algorithm.  To 
avoid  the  energy  increase  caused  by  the  repetitive  application  of  reflection  calculation,  we  extend  the 
model  in  vertical  direction  from  free  surface,  and  the  extended  part  has  the  parameters  of  background 
medium.  The  complex  screen  method  will  be  applied  to  such  an  extended  model  for  elastic  wave 
simulation.  For  each  forward  step  we  apply  the  reflection  calculation  to  the  incident  field  (only  to  upgoing 
waves)  to  get  the  reflections  and  conversions  from  free  surface.  Figure  1  also  shows  such  a  process  by  an 
upgoing  ray  from  source.  Knowing  the  reflections  and 


Figure  1.  Illustration  of  the  screen  method 


conversions,  we  can  calculate  the  scattered  field  of  the  reflected  fields  using  the  complex  screen  method 
and  then  combine  the  scattered  field  into  the  incident  field  to  obtain  new  incident  field.  In  the  following,  a 
simple  description  of  the  complex  screen  method  and  reflection  coefficient  formulas  will  be  given. 

COMPLEX  SCREEN  APPROXIMATION  FOR  ELASTIC  WAVES 


A  complete  derivation  of  the  complex  screen  method  for  modeling  elastic  waves  can  be  seen  in  Wu's 
papers  (Wu,  1994,  1996).  A  short  review  of  the  method  for  forescattering  is  given  in  this  section.  Suppose 
that  the  parameters  of  the  elastic  medium  and  the  total  wave  field  can  be  expressed  as 


and 


P(x)  =  Pq  +  8p(x) 
a(x)  =  a0  +  <5a(x) 
/3(x)  =  p0  +  <5/?(x) 

m(x)  =  u0  +  U(x) 


where  p0,  a0  and  /3(l  are  the  parameters  of  the  background  medium,  8 p(x),  <5a(x)  and  <S/3 ( x )  are  the 
corresponding  perturbations,  u0  and  U(x)  are  the  incident  field  and  the  scattered  field,  x  =  xe  x+ze  7  is  a 
2-D  position  vector  in  Cartesian  coordinates.  We  take  et  to  be  the  primary  propagation  direction 
perpendicular  to  the  thin  slab.  The  incident  field  m0  at  x0 ,  i.e.,  the  entrance  of  the  thin  slab,  can  be 
decomposed  into  plane  P  and  S  waves 

Mo(x)  =  — r  \dk\iQ  (k,x0)  +  UQ(k,x0)\elkz  (1) 

4?r  J 


where  k  is  the  transverse  wavenumber  of  plane  waves.  Superscripts  P  and  S  denote  P  and  S  waves, 
respectively.  The  forward  propagated  field  is  composed  of  primary  wave  and  forward  scattered  P  and  S 
waves.  Therefore,  at  x1  the  exit  of  the  thin  slab,  it  can  be  expressed  as 


Uf(z,  Xx) 


where 


up{z,x i)  =  (k,  x0)  +  upp  (k,  xx )+UfP  (k,x{)\ 

uSf  (z, Xj)  =  ( k , x0)  +  iiSfS  ( k , X])  +  uPS  (k,  xt)} 


(2) 


(3) 

(4) 


2  21/2  2  21/2 

where  ya  =  (ka  — k  )  and  y p  =(kp  -k  )  are  the  x  components  (propagating  wavenumber)  of  P 
and  S  wavenumbers  in  the  background  medium.  ka  =  cola  and  kp  =  of/3  are  P  and  S  wavenumbers,  a 
and  p  are  P  and  S  velocities  of  the  background  medium.  Ax  =  x1  -  x0  is  the  thickness  of  thin  slab. 
Subscripts  /  denotes  forward  scatterings,  and  PP,  PS,  SP  and  SS  indicate  the  scattering  between  different 

pp  pc  cp  SS 

wave  types.  The  scattering  and  conversion  coefficients  U y  ,  U  f  ,  U  j  and  U  j  can  be  calculated  by  the 
complex  screen  approximation  (Wu,  1994,1996). 


REFLECTIONS  FOR  P-SV  WAVES 


Suppose  that  a  plane  P  (or  S)  wave  is  traveling  in  the  background  medium  at  a  angle  of  i  (or  j )  with 
respect  to  x  (see  Fig.  1).  The  reflection  coefficients  of  the  displacement  at  the  free  surface  can  be  written 
as  (Aki  and  Richard,  1984) 
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where  z  rree  denotes  the  free  surface  location,  p  is  ray  parameter  which  is  ^2  for  P  wave  incidence,  and 

«  o 

Sln  ^  for  S  wave  incidence.  The  symbol  “  a  ”  and  “  v  ”  denote  upgoing  and  downgoing  waves, 

Po 

respectively.  From  the  screen  method,  the  incident  P  wave  at  x  =  x0  can  be  expressed  as  a  superposition  of 
plane  waves 


Uq  (z,  x0 )  =  -2—  f  dkiiQ  (, k ,  xQ)e'kz  (9) 

4  7T~  J 

Applying  the  reflection  coefficient  formulation  (5-8)  to  the  upgoing  waves  in  eq.  (9),  the  reflected  P  and  S 
waves  can  be  expressed  by 


m/(z,x0)  =  -2—  frft  1  up(k,x0 )  1  PP^-e  lkz 

4kz  J  ka 

(10) 

upp(z,x 0)  = - 2r-  f dk  1  Uq  (&,x0)  1  PP^—e~lkz 

4 J  ka 

(11) 

1  (*  d  A  V 

zzf5(z,x0)  = — y  JA:  1  Wq  (A:,^o)  1  P5 - e~lk  z 

4 n~  J  kp 

(12) 

uPS(z,x0)  =  -2-  |\flt  1 1 /(fc,x0)  1  PS^-e~ikPSz 

4 n  j  kp 

(13) 

PS  PS  12  2  2 

where  k  is  transverse  wavenumber  of  the  converted  S  wave,  using  Snell's  law,  k  =  ~k^+k  . 

I  Mq  (k,x 0)  I  is  the  amplitude  of  incident  P  wave  with  transverse  wavenumber  k  .  In  the  same  way,  the 
reflected  P  and  S  waves  due  to  the  incident  plane  S  wave  can  be  obtained  by 


Us/(z,x0)  =  -‘-y  \dk  1  usQ(k.x0)  1  SfV1"2 

4k  J  ka 

(14) 

cd  l  r  c  A  v  k$p  i  sp 

U?p(z,x 0)  = - T  \dk\uQ  (k,x0)  1  S  P - e  lk  z 

4k~  J  ka 

(15) 

u^s(z,x0)  =  — \dk  1  uq  (k,x 0)  1  S S-^—e~lkz 

4 n  J  kp 

(16) 

u?s(z,x0)  =  -L-  fdk  1  Uq  (k,x0)  1  SS^e~ikz 

4k  ~  J  kp 

(17) 

op  op  12  2  2  S 

where  k  is  the  transverse  wavenumber  of  the  converted  P  wave,  k  =  Jka  -kp  +k  .  \  u0  (k,x0)\  is 

the  amplitude  of  incident  S  wave  with  transverse  wavenumber  k  .  We  can  calculate  all  of  the  reflected 

p  s 

waves  using  eqs.  (10-17),  once  the  incident  fields  u0  and  u0  are  known.  It  is  not  difficult  to  see  that  eqs 
(10-11,  16-17),  i.e.,  the  reflected  waves  of  the  same  type  waves,  can  be  implemented  by  FFT.  However,  eqs 
(12-15),  i.e.,  the  reflected  converted  can  not  be  directly  implemented  by  FFT  because  the  nonlinear 

relationship  exists  between  k  and  kps  for  P-S  (or  ksp  for  S-P).  Although  we  can  obtain  uniform  samples 
with  respect  to  k  and  k  ps  (or  k  SF )  by  complex  variable  interpolation,  the  noise  due  to  the  interpolation  is 
so  strong  that  the  accumulated  errors  increase  fast  for  multiple  step  propagation.  In  this  study,  the  direct 
calculations  of  the  summation  integrals  are  performed  for  the  converted  reflections.  The  implementation  of 
eqs  (12-15)  by  summation  is  not  much  slower  than  by  FFT  because  only  limited  summation  is  required 
where  the  conversion  coefficients  are  not  zero. 


NUMERICAL  EXAMPLES 


To  test  the  accuracy  and  ability  of  the  complex  screen  method  used  as  a  propagator  for  crustal  wave  guides 
with  flat  surfaces,  several  numerical  examples  have  been  conducted  and  compared  with  those  calculated  by 
wavenumber  integration  method  and  finite  difference  method.  For  the  screen  method,  the  spatial  intervals 
in  x  and  z  directions  are  Ax  =  A z  =  0.25 km .  To  reduce  the  effects  of  boundary  reflection  and  large-angle 
waves,  two  smooth  window  functions  are  applied  in  space  domain  and  wavenumber  domain,  respectively. 
The  vertical  components  of  displacement  in  all  figures  shown  in  this  paper  are  multiplied  by  2. 

Figure  2  shows  synthetic  seismograms  calculated  by  the  screen  method  (solid)  and  wavenumber  integration 
method  (dashed)  for  an  elastic  half-space,  (a)  is  the  vertical  component  of  the  displacement  and  (b)  is  the 
horizontal  component.  The  homogeneous  elastic  half-space  has  a  =  4  kml  s ,  (1  =  2.5  km  Is  and 

p  =2.5 g  /  cm  .  An  explosive  source  is  located  at  16km  depth  and  has  dominant  frequency  of  1Hz.  The 
size  of  the  extended  model  used  is  1024x528  .  The  first  4  receivers  are  placed  along  free  surface  and  have 
the  offset  of  100 -124km,  and  the  last  5  receivers  are  placed  in  vertical  direction  and  have  the  vertical 
depth  0  -  32 km  .  From  Figure  2  the  reflection  and  conversion  by  free  surface  are  in  excellent  agreement.  It 
is  also  shown  that  the  use  of  direct  summation  over  incident  waves  for  converted  waves  (P-S  or  S-P)  avoids 
the  strong  noise  caused  due  to  the  complex  variable  interpolation  in  wavenumber  domain. 

Figure  3  shows  synthetic  seismograms  along  a  vertical  profile  at  the  distance  of  80 km  calculated  by  the 
screen  method  (solid)  and  wavenumber  integration  method  (dashed)  for  a  2-layered  half-space  model 
(56 km).  The  parameters  for  two  elastic  layers  are  al=6.8km/s,  ^=3.5 km/s,  pl=3.5g/cm 3, 


a2  =8.16  kml  s  ,  (32  —4.2bn/s  and  p2  =3.5  g  /  cm3 .  An  explosive  source  is  28km  depth  and  has  a 
dominant  frequency  of  0.5 Hz  .  The  size  of  the  extended  model  is  1024x320 .  From  figure  3  the  results 
from  the  two  methods  are  in  excellent  agreement  up  to  wide-angle  reflections.  The  reflections  from  free 
surface  are  angle-dependent  and  agree  exactly  with  those  calculated  by  wavenumber  integration  method. 

The  multiple  reflections  of  S  waves  behind  travel  time  of  40 sec.  from  interface  have  near  90°  incident 
angle  to  screens  and  can  not  be  modeled  accurately  by  the  screen  method.  However,  the  multiples  can  not 
survive  in  the  wave  guide  and  do  not  contribute  much  to  Lg  waves. 
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Figure  2.  Synthetic  seismograms  calculated  by  the  screen  method  (solid)  and  wavenumber 
integration  (dashed)  for  an  elastic  half-space,  (a)  is  the  vertical  component  of  the  displacement,  (b)  is 
the  horizontal  component.  An  explosive  source  is  located  at  16 km  depth  and  has  a  dominant 
frequency  of  1  Hz  .  The  first  4  receivers  are  placed  along  free  surface,  and  the  last  5  receivers  are 
placed  in  vertical  direction. 
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FIGURE  3.  Synthetic  seismograms  along  a  vertical  profile  at  the  distance  of  80 km  calculated  by  the 
screen  method  (solid)  and  wavenumber  integration  method  (dashed)  for  a  2  layered  model  ( 56 km ). 
The  source  depth  is  28 km  and  its  dominant  frequency  is  0.5 Hz  . 


Figure  4  shows  a  laterally  varying  crustal  model.  The  geometry  of  the  model  and  the  parameters  are  also 
given  in  Figure  4.  The  source  is  located  at  16 km  depth  and  has  a  dominant  frequency  of  1  Hz  .  The 
comparison  of  synthetic  seismograms  along  vertical  profile  at  the  distance  of  200 km  calculated  by  the 
complex  screen  method  (solid)  and  finite  difference  method  (dashed)  is  shown  in  Figure  5.  Excellent 


agreement  can  be  seen.  While  finite  difference  method  took  51  hours  on  a  single  processor  of  SUN 
Workstation,  the  complex  screen  method  took  only  0.6  hour. 
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Figure  4.  A  thinning  crustal  model 


Vertical  displacement  x= 200km 
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Figure  5.  Synthetic  seismograms  along  a  vertical  profile  at  the  distance  of  200 km  calculated  by  the 
screen  method  (solid)  and  finite  difference  method  (dashed)  for  a  thinning  crustal  model  ( 32 km ).  The 
source  depth  is  16km  and  its  dominant  frequency  is  0.5 Hz  . 


CONCLUSIONS 


This  paper  is  part  of  the  study  aimed  at  development  and  application  of  a  new  wave  propagation  and 
modeling  method  for  regional  waves  in  heterogeneous  crustal  wave  guides  using  one-way  wave 
approximation.  In  this  paper  we  extend  the  complex  screen  method  to  2D  P-SV  problems  for  regional 
phase  simulation  by  incorporating  the  reflection  coefficient  formulas  on  free  surface  into  the  method.  The 
comparison  of  synthetic  seismograms  calculated  by  the  complex  screen  method  and  by  wavenumber 
integration  method  and  finite  difference  method  shows  feasibility  and  validity  of  this  extension.  For  a  test 
model  with  a  propagation  distance  of  80 km  and  a  dominant  frequency  of  0.5 Hz  ,  finite  difference  method 
took  51  hours,  while  the  complex  screen  method  took  0.6  hour.  The  complex  screen  method  is  85  times 
faster  than  finite  difference  method.  For  longer  propagation  distance  and  higher  frequencies,  the  factor  of 
saving  could  be  huge.  By  this  extension,  our  ability  of  simulating  regional  phases  will  be  greatly  enhanced. 
It  makes  us  possible  to  simulate  path  effects  in  different  regions  for  various  discriminants  in  the  monitoring 
system,  such  as  Pn  /  Lg  ,  Sn  / Lg  ,  etc. 

Key  Words:  seismic  wave  propagation,  crustal  structure,  discrimination 
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